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Abstract 

In the computation of discontinuous solutions of hyperbolic conservation laws, 
TVD (total-variation-diminishing), TVB (total-variation-bounded) and the recently 
developed ENO (essentially non-oscillatory) schemes have proven to be very useful. 
In this paper two improvements are discussed: a simple TVD Runge-Kutta type 
time discretization, and an ENO construction procedure based on fluxes rather 
than on cell averages. These improvements simplify the schemes considerably - 
especially for multi-dimensional problems or problems with forcing terms. Prelimi- 
nary numerical results are also given. 
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I. introduction. In this paper we are interested in solving the system of 
hyperbolic conservation laws 
d 

(1.1a) + ^2 /,(«)*, =0 (or = g{ u, x, t), a forcing term) 

»'=i 

(1.1b) «(x,0) = tto(z) 

Here « = (tij , . . . , u m ) r , x = (xj, xj, . . . , Xd), and any real combination of the Jaco- 

d 

bian matrices ^ has m real eigenvalues and a complete set of eigenvectors. 

»=i 

On a computational grid Xj = j • Ax, t n — nAt, we use u* to denote the 
computed approximation to the exact solution u(xy, t n ) of (1.1). 

We also use the abstract form 


(1.2) ut = £(v) 

in place of (1.1a). Here C is a spatial operator. 

As is well known, the solution to (1.1) may develop discontinuities (shocks, con- 
tact discontinuities, etc.) even if the initial condition «o(^) in (1.1b) is a smooth 
function. Traditional finite difference methods, even if linearly stable, often give 
poor results in the presence of shocks and other discontinuities. Recently there 
has been a lot of activity geared towards constructing efficient finite difference 
approximations to (1.1). These include TVD (total-variation-diminishing), TVB 
(total-variation-bounded) and ENO (essentially non-oscillatory) methods. See, e.g. 
[2], [3], [4], [5], [6], [9], [10], [12], [13], [14], and the references listed therein. Many 
of the ideas can be traced back to Van Leer’s work in [15], [16]. 

Usually, rigorous analysis (e.g. total- variation stability, convergence) is only 
done for the scalar, one-dimensional nonlinear case (i.e. d = m = 1 in (1.1)). Some 
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partial theory (e.g. convergence for first order monotone schemes, and maximum 
norm stability for higher order TVD schemes) exists for scalar multi-dimensional 
problems (d > 1 in (1.1)), but a full convergence theory for multidimensional non- 
linear systems appears to be extremely difficult. However, numerical experiments 
for multi-dimensional problems and/or for systems of equations, using direct gener- 
alizations of TVD, TVB and ENO schemes give very good results. Again, see, e.g., 
[2], [4], [6], [10], [11]. We now shall confine our discussion at first to this one space 
dimension, scalar case. Systems and multi-dimensional problems are discussed at 
the end of Section 3. 

We shall always use conservative schemes of the form 

(1.3) *<4 + * a=£ 

with a consistent numerical flux 

(1.4) /y+i tii+fc); /(«,...,«)= /(«) 

in order to guarantee that any convergent bounded a.e. subsequence has as its limit 
a weak solution of (1.1), (Lax-WendrofF Theorem [8]), i.e. we construct so-called 
“shock capturing methods” . 

The total variation of a discrete scalar solution is usually defined by 

(1.5) IV(«) = 5> y+1 -« y | 


We say the scheme is TVD if 


( 1 . 6 ) 


TV{u n+1 ) < TV (u n ) 


-3- 


and TVB in 0 < t < T if 

(1.7) TV(u n ) < B 

for some fixed B depending only on TV(u°), and for all n and At such that 0 < 
nAt < T. 

A nice theoretical advantage of all TVD or TVB schemes is that they have con- 
vergent subsequences as Ax — ► 0, and, if a further “entropy condition” is satisfied, 
then they are convergent. See, e.g. [3]. 

The formal “order of accuracy” in this paper is in the sense of local truncation 
errors, i.e. if local truncation error is 0(Ax r+1 ) in smooth regions, we say the 
scheme is (formally) r-th order accurate. See, e.g. [2]. 

There are many TVD schemes constructed in the literature (e.g. (2], [3], [9], 
[10], [14]). In [10], TVD schemes of very high spatial order (up to 15th order) 
were constructed. These schemes can be used for steady state calculations (e.g. 
implemented with the TVD Runge-Kutta type time discretizations with large CFL 
numbers in [13]) or for time dependent problems, equipped with a multi-level TVD 
high order time discretization in [13] or with a Runge-Kutta type TVD high order 
time discretization in Section 2 of this paper. These are perhaps the highest order 
TVD schemes existing at present. However the definition of total variation (1.5) 
implies that these methods must degenerate to first order accuracy at extremacy. A 
TVB modification of such schemes which recovers global high order accuracy even 
at critical points is obtained in [12]. 

The above mentioned TVD and TVB schemes use a fixed, wide stencil (for the 
15th order scheme, the stencil is 17 points wide), thus restricting the- advantage of 
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going to higher order through smearing of discontinuities and resulting degradation 
of the accuracy. Numerically we observed that third order schemes work quite well 
[12], but we lost accuracy in a fairly large region near discontinuities by using a fifth 
order method. Recently Harten, Osher, Engquist and Chakravarthy constructed 
ENO schemes which are of globally high order accuracy in smooth regions and 
which use adaptive stencils, thus obtaining information from regions of smoothness 
if discontinuities are present. These methods achieve high order accuracy right up 
to discontinuities. Analysis and numerical experiments are found in [5], [6], [4]. 
At present, a convergence theory (e.g. TV boundedness) for ENO schemes is still 
unavailable. 

There are two natural directions in which to simplify the ENO or TVD, TVB 
schemes, especially for multi-dimensional problems or problems with forcing terms: 

(1) Time discretization. Usually, semi-discrete (method of lines) versions of 
ENO or TVD, TVB schemes are much simpler than the fully discrete ones. There 
are then mainly two ways to discretize in time. One is of Lax-Wendroff type, i.e., 
by using u t = -f x , u tt = (/'/*)*, •••, «y +1 = «y + Af(u t )? + ^-(u«)? + • • •, 
and then by discretizing the spatial derivatives. Many second order TVD schemes 
(e.g. Harten’s in [2]), and the ENO schemes in [5], [6], [4], used this type of 
time discretization. The main disadvantages to the procedure is that it is com- 
plicated to program, especially for multi-dimensional problems with forcing terms. 
One can see this by writing out a third order approximation to the equation 
u t + /i(«)*i + /2(*0*, = g{u, xi,xa, t). Moreover it is not easy to prove that 
this results in a TVD or TVB method, even if the original method-of-lines ODE 
and its Euler forward version are both TVD or TVB. The numerical results have 
proven satisfactory, but, speaking theoretically, only second order in time TVD or 
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TVB schemes exist of this type. Another way to discretize in time is to use a multi- 
level or Runge-Kutta type ODE solver. This is much simpler to program than the 
Lax-WendrofF type of discretization for multi-dimensional problems or for problems 
with forcing terms, so it is widely used for numerically implementing a method of 
lines approximation. However, usually only linear stability analysis is available in 
the literature, which is certainly not enough for our purpose since linear stabil- 
ity does not imply convergence if shocks or other discontinuities are present. This 
is particularly true for ENO schemes which use moving stencils. Linear stability 
analysis is based on the fact that the stencil is fixed and the error accumulates in 
a predictable pattern, hence it does not apply to ENO schemes at all. For these 
reasons we consider TVD time discretizations . In [13], a class of multi-level TVD 
time discretizations were constructed and analysed, (numerical results can be found 
in [12]). However, for easy starting and for storage considerations, one step Runge- 
Kutta type schemes are preferable to multi-level methods. In Section II of this 
paper we present a class of high order TVD Runge-Kutta type time discretizations. 

(2). Avoiding the using of cell-averages. The ENO schemes constructed in [5], 
[6], [4] are for cell- averages but involve point values as well. Hence a reconstruction 
procedure is needed to recover point values from cell averages to the correct order, 
which can be rather complicated, especially in multi-dimensional problems. It is 
desirable to use the moving-stencil idea directly on fluxes to get ENO schemes 
without using cell-averages. In Section III of this paper a class of such ENO schemes 
is constructed. 

Some encouraging numerical results obtained by using schemes constructed in 
this paper are included in Section IV. 


- 6 - 


We conclude these introductory remarks by noting that R. Sanders [17] has 
recently devised third order accurate TVD methods which degenerate to second 
order at extrema. He defines the variation of the numerical solution as the variation 
of an appropriately chosen piecewise parabolic interpolant. The numerical results 
are very good. However this technique has no method of lines analogue, so we omit 
it from our present discussion. 

II. High Order Runge-Kutta Type TVD Time Discretizations. Define 

(2.1) w = T(u) = (I + AtL){u) 

where T and L are nonlinear discrete operators, L is a r-th order discrete approxi- 
mate to the spatial operator C in (1.2): 

(2.2) L(u) = C{u) + 0{ Ax r ) 
if u is smooth. 

Our goal is to get a fully r-th order approximation to the differential equation 

(1.2) of the form 

(2.3) u n+1 = S{u n ) 

(The operator S depends on T). This means that if u(x,t) is an exact smooth 
solution of (1.2), then 

(2.4) u{xj, t n+1 ) - S{u n )j = 0{ Ax r+1 ). 

We also want the scheme to be TVD: 

(2.5) 7V(5(tt)) < TV{T{u)) 

under suitable restrictions on At (or, equivalently, on the CFL number A). 
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We call a time discretization (2.3) r-th order TVD if it satisfies (2.4) and (2.5). 
If the spatial operator T in (2.1) is TVD or TVB: 

(2.6a) 7V(T(ti)) < TV(u) 

or 

(2.6b) TV(T{u)) < TV(u) + MAf 

for 0 < M uniformly bounded as At — ► 0, then the fully discrete high order scheme 
(2.3) is TVD or TVB, owing to (2.5). 

In [13], a class of multi-level type high order TVD time discretizations was 
constructed. Numerical experiments in [12] were very promising . But there are two 
disadvantages of multi-level type methods: (i) for an m-th level method the first m — 
1 levels have to be calculated by other methods to the same order of accuracy (e.g. 
by using Taylor series expansions); (ii) we have to store all m level datas, creating a 
rather large storage requirement, stretching up to and beyond the limits of present 
day computers for physical problems arising e.g. in computational aeronautics. At 
present, Runge-Kutta type methods are more often used in discretizing the method- 
of-lines than are multi-level methods. Since the former consists of one-level methods, 
they are self-starting and reduce storage requirements significantly. In the following 
we will analyze the nonlinear stability (TVD) of a class of such methods. 

Assume (2.1) is TVD (or TVB) under a suitable CFL restriction 
(2.7) A < A 0 


We may also need an approximation to — L to the spatial operator which 


we take to satisfy 

(2.8) 

(2.9) 


6 = T(v) = (I - AtL)(u) 
L(u) = £(u)+0( 4i r ) 
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where (2.8) is TVD (or TVB) under the same CFL restriction (2.7). 


As an example, a very simple first order (r = 1) TVD approximation to 

u t =u x = £(u) 


is obtained via simple upwind differencing: 


L(u) = 


u(x + Ax) — «(j) 

Ax 


This scheme (2.1) is TVD if A = < A 0 < 1. 


The approximation L(u) is defined as 

r,..x _ u(z)-u(z-Az) 
L{U) ~ Ai 

and (2.8) is TVD again for A = ^ < A 0 < 1. 


This procedure (2.8), (2.9) easily generalizes to any conventional TVD, TVB, 
or ENO- approximation (2.1) satisfying (2.2). 

The general explicit Runge-Kutta method (we use explicit methods to avoid 
solving nonlinear equations) for (2.1) is 


(2.10a) 

«-l 

y(») =|i (0) _|_ A/ ^ CilfL(U ^), 

M , 

II 

to 

3 


fc=0 


(2.10b) 

u(°) = u n , « (m) = u n+l 



If the operator L also depends explicitly on t, as is the case when the forcing 
term g in (1.1a) depends explicitly on t, or when we have time-dependent boundary 

conditions, the general explicit Runge-Kutta method takes a more complicated form 

*’— 1 

(2.11a) «<*> = u + At ^ CikL{u^ k \ + d^At) 

fc =0 


where 
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fc-1 

(2.11b) d k = J2 c “ 

1=0 

For details, see any numerical ODE text, e.g [1], [7] (any such method is usu- 
ally written in a slightly different form, using hi, k< 2 , . . • , but that form is clearly 
equivalent to (2.10) or (2.11)). 

We shall restrict our attention to (2.10); the generalization to (2.11) is clearly 
straightforward, using (2.11b). 

In order to get conditions for TVD, we rewrite (2.10) as follows: For a,-* > 
»— l 

0, = 1, we have 

k= o 

= ^2 a ik u(0) + At ^2 c ikL{u (k) ) 

k= 0 fc=0 

1 k— 1 

= a t0 « (0) + ^2 a «fc( u(<) “ c kiL{uW))+ 

k=l i=o 

+ At^Tc ik L{u( k) ) 

k- 0 

= + (c,-fc - c <fc a,/)AfI(ti (fc) )] 

fc=o e=k+i 

i- 1 

So if we let /?,•* = c ik - ^ c ik aa, then (2.10) may be written in the following 

t=k+ 1 

equivalent form 

(2.12) «'•» = £>»«<*> + A»At£(«(*))] 

k= 0 


It is well-known that we can get (m + l)-th order accurate methods in the form 
(2.10) or (2.11) for m < 3; m-th order methods for m = 4, 5, 6; or (m — l)-th order 
methods for m = 7,8 (see, e.g. [1]). 
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For the classical 4-th order Runge-Kutta methods, the constants c,* in (2.10) 
are all non-negative. However, /?,•* in (2.12) may well be negative. In order to obtain 
TVD we apply a trick used in [13], i.e. we replace L in (2.12) by L in (2.8)-(2.9) 
whenever /?,-* is negative. Now (2.12) becomes a convex combination of TVD (or 
TVB) operators under the CFL restriction 

OK • jl 

(2.13) A < Aominrr— r 

v ' <.* \Pik\ 

and we easily get the following 

PROPOSITION 2.1. Scheme (2.12) is TVD under the CFL restriction (2.13), 
if L is replaced by L when /?,•* is negative. 

REMARK 2.1. The previous proposition may be put into a more general 
framework as follows. The TV in (2.5) and (2.6) (a, b) may be replaced by G(u), 
any convex mapping into the non-negative real line, where u, S(u), T(u) belong 
to a Banach space of functions B Ax . Also if u is a smooth solution of 1.2, then 
u(x,t n+1 ) - S(u(x, t n )) = 0((Ax) r+1 ). The statement in the proposition can then 
be replaced by: G(u m ) < C7(u 0 ); moreover the formal order of accuracy is still r. 

Now our goal is to choose the a,* and /?,-* such that (2.12) is of the highest 
possible order and such that the CFL restriction (2.13) is optimal. We would also 
like to minimize the number of negative /?,-*’ s in order to reduce the computational 
work involving L. 

One easy way to do this is to use a standard Runge-Kutta method and then 
rewrite it in the form (2.12) to get a,* and flu-- Unfortunately, most classical 
Runge-Kutta methods lead to small CFL numbers in (2.13) as well as negative 
Pik s. Hence the best way is to consider (2.12) directly. Straightforward but tedious 
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Taylor expansions and an analysis of possible parameters (which we omit) leads us 
to the following results: 

i) Second order case, m = 1 . 


For accuracy 


(2.14) 


' £*20 = 1 ~ a 2 l 
i 020 — 1 _ 2^ ~ a 2l/?10 
^ 021 — 2/Jio 


Pio, a 2 i are free parameters. 


It can be verified that the “optimal” scheme (considering CFL restriction (2.13) 
and whether L appears) is 

{ «(!) = «( °) +AtL{u(°)) 
u (2) _ i w (o) + I^£( u ( i)> . 

CFL# = 1 

Notice that L does not appear in (2.15). This is equivalent to 

( ul 1 ) = « 1 °) + A tL(u^) 

(2 ' 16) I u( 2 ) = + |A<L(u(°)) 4- %AtL(uW) 

which is the classical Heun’s method or modified Euler method [ 1 ]. 


ii) Third order case, m — 2. 


For accuracy 


(2.17) 


( a _ 3ftio—2 

^32 - 6P(/9 10 -P) 

/^21 6 /? io 0*3 

< a *— <*»3 0io 03i — P 013 

/?31 - ^ 

030 = 1 — <>31 010 ~ O32 P ~ 031 — 032 
, 020 —P- «21 PlO ~ 021 


021 , 030 , a 31 , /? 10 and P = /? 20 + £*21 /?io + A 21 are free parameters. The solution 
(2.17) is written in convenient inductive form. 
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Extensive searching leads to the following preferred scheme: 

'«(*)= «(°) + 

«( 2 ) = + ^AtL(u^) 

(2 ' 18) | «(») = l u (») + §«<»> + f AfL(u< J >) 

. CFL # = 1 

Notice that L does not appear in (2.18). Also in computing we only need 
L(u^*~^), so there is no need to store the previous £(«(*)), lowering the storage 
requirement significantly. 

(2.18) is equivalent to 

{ «(1) =U (0) +A *£( U (0)) 

= «(°) + jA^L(u^ 0 )) + ^AtL(u^) 

«( 3 ) = u(°) + gAf£(u(°)) + 5 A tL(u^) 4- |A tL{uW) 

We have been unable to identify (2.19) with any of the “classical” third order 
Runge-Kutta methods. On the other hand, the classical third order Runge-Kutta 
methods in [7], when written in equivalent forms (2.12), lead to negative and 
small CFL numbers (2.13), and are hence inferior to (2.19). 

(iii) Fourth order case, m = 3 

For accuracy we get a system of 7 equations with 16 unknowns, so there are 
9 free parameters. Unfortunately this time the solution is not easily obtained in a 
convenient form. In [1] a general solution with two parameters is given for the form 
(2.10). We can certainly rewrite it in the form (2.12). Extensive searching seems to 
indicate that we cannot avoid negative constants /?,* this time. The classical fourth 
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order Runge-Kutta method can be written in the form (2.12) as 

' = u<°) + ±A tZ(u<°>) 

«( 3 ) = |u(°) — |A tL(u^) + + ^A tL(uW) 

(2.20) i «( 3 ) = I«(°) - IA/Z(«(°)) + §«<*) - ±A tL{uW) + | uW + A tL(uW) 
«( 4 ) = |«( l ) + jAfZ/(«(^) + + ^A tL(uW) 

< CFL# = | 


Notice that we have to compute L(u (°)) and L(u^). If we use the more 
awkward definition of u ^ 


( 2 . 21 ) 


/»\ 6431 /q\ 18769 . ; , 18769 137 . i . fi\. 137 A fiK 

tt (3) = «(°) A fI(u l0, )H Af£(u ll, )H u (2 '+AfL(«i 3) ) 

80000 160000 ' '80000 400 ' ' 200 1 ’ 


then the CFL# can be raised slightly to 

Another fourth order scheme with a slightly larger CFL # is: 
u^ 1 ) = u<°) + ^AtL(u^) 

u< 2 ) = §u(°) - | A tl(u(°)) + fu* 1 * + fAfL^ 1 )) 

« (3) = 5§§Jo« (0) ~ m^i(uW) + JjgL U {i) - giA tl[nM)+ 
+g«(3) + |AfL(«< 2 )) 

ul 4 ) = | A fL(u^) + + %AtL(uW) 

CFL # = 0.87 

We still need to compute L(u (°)) and L(u^). 

(iv) Fifth order case, m = 5. 



We simply write out the form (2.12) corresponding to a fifth order method 
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given on page 143 of [7]: 


(2.23) 

' «<*) = «l°> + | A tL(uW) 

«( 2 ) = |«(°) + ^i^ 1 ) + |A tL(u^) 

u< 3 ) = |«(°) - ±At£(u(°)) + + |A tL(uW) 

«( 4 ) = |u(°) - fjAtL(u(°)) + — gf A tL(uW) + + |A tL(uW)+ 

< £u< 3 ) + £Af£(u< 3 )) 


m, 

12000 ( 


„(•) = }„(«) + - Aa (£(«<*>) + A„(3) + |AiI(«' 3 ') + if«< 5 > + A A tL(«W) 

( CFL# = A 


Notice that we need to compute L(u(°)), L(u 0)), Z(«( 3 )). 


III. A Simplified Version of ENO Schemes. To use the Runge-Kutta 
type TVD time discretizations in Section II, we must have a spatial discrete opera- 
tor (2.1) to start with. Theoretically one would like to use a TVD or TVB operator 
T satisfying (2.6), because then the full scheme (2.3) would be TVD or TVB. But 
as indicated in section I, the existing high order TVD or TVB schemes may smear 
discontinuities and pollute the solution (i.e. we may not get high order accuracy in 
a fairly large region near discontinuities), due to the fixed, wide stencil. The ENO 
schemes constructed in [5], [6], [4] are very promising experimentally and appealing 
conceptually, but the fact that they use cell-averages as well as point values via a re- 
construction procedure, and that they were implemented using a Lax-Wendroff type 
time discretization, makes them rather complicated to program, especially in multi- 
dimensional problems, or problems with forcing terms. The Runge-Kutta type TVD 
time discretizations in Section II equipped with semi-discrete ENO schemes will sim- 
plify them in many cases (although there is no rigorous theory concerning TVB of 
semi-discrete ENO schemes or their Euler forward version, analysis in many cases 
and numerical experiments strongly support that the total variation increase at each 
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step is 0(Ax r ) for a r-th order ENO scheme, see [6]. Hence the full scheme (2.3) 
in this case should also be TVB). In this section we further simplify non-oscillatory 
methods by deriving a version of ENO schemes using only fluxes, not cell-averages. 

We start with a simple first order monotone Lax-Friedrichs type of scheme. If 
we define 

(3.1) /+(«) = i(/(«) + an), /"(«) = i(/(u) - au) 
where a > max|/'(«)| is a constant, then clearly 

(3.2) f +, (u) > 0, /“'(«) < 0 

(3.3) /*M+r («)=/(«) 

The Lax-Friedrichs scheme is simply (1.3) with the numerical flux defined by 

(M / J+} =//+/,- +1 


Taylor expansion reveals the existence of constants a 2 , a 4 , . . . , a 2 m - 2 , • • • , such 
that if 


m— 1 


(3.5) 


0* 


f j+i = f i+i + £ ^ k (g-sf) i+i + 0(A^ m+I ) 




then the scheme (1.3) will be 2m-th order accurate in space in the sense of (2.2). 


For example, a 3 = — 04 =-^,.... 


In light of (3.4), it is natural to require 


(36) /, +} = /; +i +/- +i 

and to define the positive flux ft x and the negative flux h . (in the meaning of 

J'T' 2 J T 2 

(3.2)) separately. 
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For accuracy, we require ff +x and f J+x to satisfy (3.5) separately: 


•S+t 

m — 1 


(3.7) 


glk 


fU = /Ai + E 


k= 1 


We achieve (3.7) by using polynomial interpolants p * of /* to the correct 
order: 

(3.8) p*(x) = /*(«(*)) + 0(A* J ">+‘) 

near x = x J+ x, then define 


m— 1 


(3.9) 


5“ 


/,*+} = Py (*/+)) + £ <*»( A *)“(flp*l>/)._. 

k= 1 ^ 


Clearly if (3.8) is true, then the fluxes /£ A defined by (3.9) will satisfy (3.7). 

It is in constructing the interpolating polynomials pf{x) that we use the ENO 
moving stencil idea: in order to achieve (3.8) pf(x) can be polynomials of degree 
2m interpolating /*(«(!:)) at any 2m + 1 points near x J+ i. We use the ENO ideas 
in [5], [6] and [4] to choose the 2m + 1 points automatically from the smoothest 
possible region, but start with the correct one according to (3.4). 


The algorithm can be written as follows: For constructing (x) : 

(3.10) (1) k^ n = kSL=j, <?! h 0) (x)=/ + («y) 

(2) Inductively, assume we have Araax* 1 and then we com- 

pute the n-th divided differences of / + (u(x)): 


(3.11a) 

(3.11b) 


fl(n) = f + W{x k 

b^=f + [u(x k 


(n-l)),...,«(x.(n-l) )] 

min *n»ax * * 

(n-i) )j . . . , U(X, (n-1) )]. 

«n I n * "milt 


min 
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We proceed to add a point to the stencil according to the smaller n-th divided 
difference: 

(i) If|a (n >|>|& (n >|, then 

(3.12a) =*>(") 

(-1 iow u( n ) — _ i fc(») = 

(O.IZD) K min *min A > ft max ^max 

(ii) If |a(")| < |6( n >|, then 

(3.13a) c (n) =a( n ) 

(3.13b) *<& = *£" l) , + 1 

and finally 

(3.14) Q^(x) = Q { ?~ l) {x) + c (n) n^ ( l ) . 1 ,(x- x h ) 

min 

(3) p+( *)=<4 Sm ’(*) 

For constructing pj(ar): 

(!) = ; + 1, QL 0) (s) = /“(uy+i) 

(3.15) (2) same as (2) above with / + replaced by f~ and Q+ replaced by Q _ 
(3) P j(x) = qL 2 ""(i) 

REMARK (3.1). 


(a) For the first order scheme we just get back the Lax- Friedrichs flux (3.4); 


- 18 - 


(b) The second order scheme here is similar to the usual minmod second order 
TVD scheme, [9] except that the minmod function is replaced by choosing the value 
closer to zero (we omit the details of the derivation here); this is still a TVD method. 

(c) For the linear equation «t + au x = 0 (the scheme is still nonlinear!), the 
schemes here are equivalent to the ENO schemes in [6] using the primitive function 
reconstruction, except for a possible difference in the choice of stencil. The ENO 
schemes in [6] choose the stencil according to the divided difference table of cell 
averages of u, rather than that of / ± (u) here. We again omit the details of derivation 
here. Since the ENO schemes in [6] worked so well numerically, and our simplified 
schemes here are equivalent to those in [6] for linear equations, we expect ours to 
work as well. For preliminary numerical results see section IV. 

As mentioned before there is at present no rigorous theory about TVB of this 
type of ENO schemes. Here we make two observations for our schemes along these 
lines: 

(1) For smooth solutions all the divided differences (3.11) should be bounded 
(by the maximum norm of the n-th derivative of / ± , times some constant). So if 
we use 

(3.16a) c( n > = min(|c( n )|, M^) sign (c^) 

or 

(3.16b) c* n ) = min(|c< n )|, M< n > Ax n ~ 2 ) sign ( c ( n >) 

in the place of c in (3.14) for n > 2, where M are constants which are related 
to the maximum norm of the n-th derivative of in initial smooth regions, it does 
not affect the accuracy in regions of smoothness. We then get a TVB scheme. We 
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can easily see that our flux with (3.16) satisfies 

(3.17a) f i+ i=fZT + e i+i 

where 

(3.17b) |c y+ ^| < MAx 3 

with the constant M depending on the , and /T^ 3 is the second order TVD 
flux mentioned in Remark (3.1b) above. Now (3.17) clearly implies TVB of the 
scheme. Numerically we do not see any essential difference by using or not using 
(3.16), hence we strongly believe that ENO schemes without (3.16) are also TVB, 
at least for most practical problems. 

(2) The approach in (3.1) - (3.15) is not the only possible one which can be 
used to construct an ENO scheme based on interpolating fluxes. At an early stage 
of our current work we used another approach: starting from / + and f~ in (3.1), 
then for constructing p+ (x): 

(!) frJin = j ~ 1. *!i!* = j, <4 1} (*) = f + {Uj) + / + [u(xy_ 1 ),tt(Xj)](x - Xy); 
(2) Sc (3), same as the procedures (2) and (3) above in (3.10), 

Similarly, for constructing pj(x ): 

(1) fciun = frmic = / + 1, (^(z) = /“(«,*) + /~[«(zy), (zy+i)](z - xy); 


(2) and (3), same as the procedures (2) and (3) above in (3.15). 

Then we take py(x) = p+ (x) + pj(x), and write our scheme as 

«/ +l = - At{^pj{x)) x=x . 


(3.18) 
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Notice that the scheme (3.18) is simpler than (3.9) - (1.3). The only trouble is that 
it is not in conservation form (1.3). However, if we use (3.16), then it is easily seen 
that (3.18) can be written as 

(3.19) «y +1 = - Hfj+i ~ ff-i) + Ax 3 c y 

where |cy| < M, and /f^ A is the first order Lax- Friedrichs flux (3.4). We can call 
such schemes “essentially conservative” because the most important property of a 
conservative scheme-the conclusion of the Lax-WendrofF theorem in [8] - is still 
valid. Since the scheme deviates from a first order monotone scheme (not only a 
TVD scheme) by M Ax 3 , we have even a stronger theory than before - we have 
the entropy condition, hence full convergence (not just of a subsequence), and also 
convergence in multi-dimensional scalar problems i.e. we have every convergence 
property first order monotone schemes have. Unfortunately numerical experiments 
indicate that in some cases, (3.18) is inferior to the fully conservative (3.9) - (1.3). 
An illustrative example is to compute the Riemann problem for Burgers’ equation 
u t + uu x = 0 with a moving shock (e.g. «i e ft = ^ bright = — 1) using the fifth 
order versions of (3.18) and (3.9) - (3.15), (1.3), equipped with a fifth order multi- 
level TVD time discretization in [13]. The procedure (3.9) - (3.15), (1.3) gives good 
results, with or without (3.16); while without (3.16), the non conservative (3.18) 
gives the wrong shock location. With (3.16) the shock location becomes correct, but 
the mechanism that enforces this causes a rather severe smearing of the shock. For 
these reasons we abandoned the simple and theoretically pleasing version (3.18). 

Finally let us point out that the Lax- Friedrichs building block is only a con- 
venient one; we may also use other monotone or E-fluxes (see, e.g. [10]) as our 
building blocks. Of course it is not always possible to assoicate / + and f~ as in 
(3.1) with each E-flux such that (3.2), (3.3), (3.4) is valid, but careful inspection 
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reveals that we do not need to use the values of / + and / - only their divided 

differences. For each E-flux h J+ x we can define 

(3.20) x = fj+i ~ x , df j+i = ^y+ a — /y 

where d/t, A and d/7 x replace the first (undivided) differences ff +l — /+ and 
/r ) _ 1 — /r. Hence we can just use the divided difference tables of d/ + and d/ - in 
place of the divided difference tables of f + and f~ in constructing p+ and p ~ , and 
define /+, fj in any consistent way such that f+ +fj = /y, e.g. /+ = /y, /y* = 0. 
By (3.20) 

(3.21) df* +i +dfr +i = f i+l - fj 

hence if and pj use the same stencil then pf{x) +Pj{ x ) is a polynomial inter- 
polating f(u(x)), thus accuracy is guaranteed with (3.9) - (3.6); if the stencils are 
different, say p+ has the same stencil as pj but p+ does not, then it is easy to show 
that p+ —pf is a sum of r-th order undivided differences of dft +k (r is the order of 
the polynomials p+, p7) hence as long as these are 0(Ax r+l ) (valid if the JF-flux 
h J+ x is smooth up to order r) we still have the correct accuracy. 

The reason one might consider general F-fluxes as building blocks is that the 
Lax-Friedrichs flux is considered to be too dissipative. While the first order Lax- 
Friedrichs scheme is much inferior to upwind schemes (e.g. to Godunov’s or the 
Engquist-Osher schemes), our numerical experiments show that higher order ENO 
schemes using Lax-Friedrichs building blocks work quite well (although they are still 
slightly inferior to the same order ENO schemes based on upwind building blocks, 
the difference is much smaller than that in the first order case) . The advantage of 
the Lax-Friedrichs flux is that it is C°°, hence the ENO schemes based on it have 
full high order accuracy. On the other hand, most other E-fluxes - Godunov’s, 
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Engquist-Osher’s, Entropy condition satisfying version of Roe’s, etc. - are not 
smooth at sonic points (points at which f'(u ) = 0), hence ENO schemes based 
on them using the methodology of this paper will lose accuracy at sonic points. 
Although we may overcome this by smoothing those i?- fluxes at sonic points, in 
most cases the simple Lax- Friedrichs building block should be good enough. 

Problems in multi-dimensions are approximated by applying the procedure 
described in (3.1) - (3.15) or its generalization (3.20), (3.21) to each of the terms 
in (1.1a). The Range- Kutta methods devised in section 2 are then used, with 
CFL numbers shrunk by a factor (d)~ x . 

Systems of equations are approximated using by now familiar field-by-field 
decompositions ideas. In (3.1) we replace the scalar constant a by a constant 
matrix al, a = {<*,y}Jy =1 , where the eigenvalues of are non-negative and of 
are non-positive. 

Obviously a might be taken to be a sufficiently large positive scalar multiple of 
the identity, but this might also lead to some smearing of discontinuities associated 
with slower waves. Other more practical choices might involve freezing §£ at some 
constant state 8, diagonalizing: 


U - m 


M«) 


Arf(u) 


T-Un) 


and letting 


a = T(u) 


T-'iu) 


Xd. 


where each A,- > |A,(u)| throughout the region. To be safe, the margin of difference 
between each A,- and the maximum value of |A,(u)| has to be sufficiently large. 
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In any case, the corresponding pf (x) are obtained with the help of the left 
and right eigenvectors of (u 3 ), which we denote by ftp and r^\ i = 1, ...d. We 
interpolate ftp • /^(u) obtaining ftp • pf(x ) exactly as in (3.10) - (3.15). We then 
define 

i-l 

The fluxes M. are defined through (3.9). 

Generalizations of the type described in (3.20), (3.21) using approximate Rie- 
mann solvers for h J+ A , appropriately smoothed at sonic points, may also be ob- 
tained. 

Work is currently under way with various colleagues applying these methods 
to Euler’s equations of compressible gas dynamics in multi-space dimensions. 


IV. Preliminary Numerical Results. The numbers in this section are often 
written in exponential form, e.g. 4.2-3 means 4.2 x 10 -3 . 


EXAMPLE 1. The ENO schemes (3.1) - (3.15) in section III combined with 
the Runge-Kutta type TVD time discretizations (2.19) - (2.20) in section II are used 
to solve the nonlinear Burgers equation with periodic initial conditions: 


«< + (t). = 0 

ttfar.O) = h + i 


- 1 < x < 1 


The exact solution is smooth up to t = J, then it develops a moving shock 
which interacts with the rarefaction waves. We get the exact solution by using a 
Newton iteration. For details, see [6]. 

Since there is a sonic point, we use the smooth LF (Lax- Friedrichs) building 
block in our ENO schemes. Both 3-3-LF-ENO (third order in time and space ENO 
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schemes with Lax- Friedrichs building blocks) and 4-4-LF-ENO are used. 

We use a CFL number of 0.8 for 3-3-LF-ENO and 0.6 for 4-4-LF-ENO. 

The errors of the numerical solutions at t = 0.3 are listed in Table 1. Since the 
exact solution is still smooth, we get the full order of accuracy in both L\ and Loo 
norms. 

At t = the shock begins to form. We use Ax = 55 and print out the errors 
at 10 points near the shock: 

3- 3-LF-ENO: -8.7-4, -2.6-3, -7.1-3, -1.3-2, -1.2-1, * 

8.2-2, 8.0-3, 1.1-3, 2.6-4, 1.9-4 

4- 4-LF-ENO: -1.5-4, -3.2-4, -1.6-3, -1.5-2, -1.2-1, * 

7.8-2, 7.9-3, 9.2-4, 1.8-4, 9.5-5 

where the * is the position of the shock. 

We see that there is a very good shock transition. (No oscillations are observed). 
Figures 1-6 show the shock transitions. 

In smooth regions the numerical solutions are very accurate. We compute the 
Li and L «> norms in the region a distance of 0.1 from the shock (i.e. |x— shock location| > 
0.1) and list them in Table 2. From the table we can see that the errors are of the 
same magnitude as in the smooth case when t = 0.3. 

At t — 1.1, the reaction between the shock and the rarefaction waves is over. 

The solution becomes monotone between the shocks. We again print out the errors 
at 10 points near the shock for Ax = 

3-3-LF-ENO: -1.0-4, 4.6-4, 4.2-4, -3.3-2, -8.3-3, * 


- 25 - 


3.9-2, 1.7-3, -1.6-4, -2,5-5, -7.0-6. 

4-4-LF-ENO: 1.1-5, 6.6-4, -1.6-3, -5.7-2, -1.3-3, * 
5.7-2, 2.7-3, -2.4-4, -5.8-5, -1.2-6 


Figures 7-12 show the shock transitions. 

The errors where the solution is smooth are again listed in Table 2. 
We can see the excellent behavior of ENO schemes in this example. 
EXAMPLE 2. A two-dimensional version of Example 1 
(4.2) < ; -2<x, y<2 


f «• + (¥). + ttf), = 0 

l «(*, y, o) = i + a sin 


is tested using the same schemes as in Example 1 in a dimension by dimension 
fashion (i.e. T(u) = (I + AtL x + A tL y )(u) in (2.1), together with the R — K time 
discretization: (2.19)-(2.20)). The exact solution is one-dimensional depending only 
on f = x + y, however our grid points are rectangular in (x, y) coordinates, and 
thus this example is a truly 2-dimensional test problem. 

The CFL number is always taken to be half of the one dimensional analog, i.e. 
0.4 for the 3-3-LF-ENO and 0.3 for the 4-4-LF-ENO. 


As in Example 1, we collect the L\ and L<» errors at t — 0.3 (smooth solution) 
in Table 3 and the Li and L<x> errors in regions at a distance of 0.1 from the shock 
at times t = % and t = 1.1 in Table 4. We also print out 10 points near the shock 
when x = 0, t = ^ and t = 1.1 for Ax = A y= 

t ~ * ’ x = 0: 

3-3-LF-ENO: -9.7-4, -2.3-3, -7.6-3, -4.5-3, -1.2-1, * 
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8.2-2, 7.9-3, 1.1-3, 2.6-4, 1.9-4 

4-4-LF-ENO: -1.5-4, -3.2-4, -1.6-3, -1.5-2, -1.2-1, * 

7.8- 2, 7.9-3, 9.2-4, 1.8-4, 9.5-5 

t - 1.1, x = 0: 

3- 3-LF-ENO: -1.0-4, 4.6-4, 4.2-4, -3.3-2, -8.3-3, * 

3.9- 2, 1.7-3, -1.6-4, -2.5-5, -7.0-6 

4- 4-LF-ENO: 1.1-5, 6.6-4, -1.6-3, -5.7-2, -1.3-3, * 
5.7-2, 2.7-3, -2.4-4, -5.8-5, -1.2-6 


The shock transition graphs are very similar to Figures 1-12, hence we omit 
them. 


We observe essentially the same results as in the 1-dim Example 1. This indi- 
cates that our ENO schemes work well in multi-dimensional problems. 

EXAMPLE 3. We use the same schemes as in Example 2 above to solve a 
linear problem 


(4.3) 


{ «t + U* + Uy = 0, — 1 < X, Jf < 1 


where S = {(x, y) : |x — y\ < |x + y| < is a unit square centered at the 
origin and rotated by an angle of * (see [4]). We use Ax = and run the scheme 
up to t = 16 (8 periods in time), in order to study the stability and the amount of 
smearing of discontinuities of these methods. 

The numerical solutions at t — 2 (after 1 period in time) and at t = 16; y — 0 


and y = -0.4, are displayed in Figures 13-20. 
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Observations: (1) the 2 dimensional schemes are stable under CFL numbers 
one half of those used for 1 dimension; (further experiments using 1-dimensional 
CFL numbers led to instability-overflow), 

(2) The 4-th order scheme resolves the discontinuities better than the 3rd order 
method. 

(3) Overshoots and undershoots, if any, are negligible. 

EXAMPLE 4. We did not prove for these ENO methods that limit solutions 
satisfy the entropy condition. However, numerical experiments in [6], including 
some tests using nonconvex fluxes, indicated the convergence of ENO schemes to 
the correct entropy solution. We test our schemes 3-3-LF-ENO and 4-4-LF-ENO 
for Riemann problems for two such fluxes. One is 

((4.4)) f(u) = i(u a - l)(u 3 - 4) 

with ul = 2, ur = —2 (the exact solution is a shock followed by a rarefaction 
wave followed by another shock) and with «£, = —3, ur = 3 (a stationary shock at 
x = 0). See [6] for details. Our schemes converge to the correct solutions in both 
cases with good resolution. The results are displayed in Figures 21-32. 

Another nonconvex flux we test is the well known Buckley-Leverett example: 

<“- 5 > '<“> = 40 +?-«)». 

with initial data u = 1 in [-£,0], and u — 0 elsewhere. 

The exact solution is a shock-rarefaction-contact discontinuity mixture. Our 
schemes resolve the correct solution well. However, the 3-3-LF-ENO (using the Lax- 
Friedrichs building block) smears more than the 3-3-EO-ENO (using the Engquist- 
Osher building block) (Figures 33-38). Although in this example f'(u) > 0 so 
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there was no need to smooth the flux, the observed improvement using the upwind 
building block indicate that sometimes it is worthwhile spending the effort to smooth 
the EO or other upwind flux rather than to use the simple LF building block. 
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Table 1. (Example 1). 

t — 0.3; E : type of error; r : numerical order 


Ol 

Loo 

L i 


3-3-ENO 

m 


B 

3-3-ENO 

r 4-4-ENO 

r 

X 

10 

2.6-3 


2.1-3 


7.2-4 

5.3-4 


X 

30 

2.3-4 

3.50 

1.1-4 

4.25 

5.7-5 

3.66 2.5-5 

4.41 

X 

40 

3.1-5 

2.89 

5.4-5 

4.35 

6.5-6 

3.13 9.2-7 

4.76 


Table 2. (Example 1). 


Errors in smooth region |x— shock| > 0.1; Ax = 



L i 

t = 

2/ir 

t = 

1.1 

t = 

2/ff 

< = 

i.i 

3-3-ENO 

4-4-ENO 

3-3-ENO 


3-3-ENO 

4-4-ENO 

3-3-ENO 

4-4-ENO 

8.7-4 

1.5-4 

1.6-4 

1 

1.6-5 

3.5-5 

8.6-6 


1.4-6 
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Table 3. (Example 2). 


t = 0.3; E : type of error; r : numerical order 


Ax = Ay\ S 

Loo 

Li 


3-3-ENO 

r 

4-4-ENO 

r 

3-3-ENO 

r 

4-4-ENO 

r 

1 

s 

2.7-3 


2.1-3 


1.4-3 


2.7-4 


_L 

10 

2.3-4 

3.55 

1.1-4 

4.21 

1.1-4 

3.67 

1.3-5 

4.40 

J- 

20 

3.2-5 

2.85 

5.6-6 

4.36 

1.3-5 

3.08 

4.5-7 

4.80 


Table 4. (Example 2). 

Errors in smooth region \(x, y) — shock) > 0.1; Ax = Ay = ^; 


Loo 

L x 

t = 2/jt 

t- 1.1 

t = 2/jt 

<=1.1 

3-3-ENO 

4-4-ENO 

3-3-ENO 

4-4-ENO 

3-3-ENO 

4-4-ENO 

3-3-ENO 

4-4-ENO 

9.9-4 

1.5-4 

1.6-4 

1.7-5 

7.9-5 

4.8-6 

1.5-5 

B 




























- 31 - 


In the following figures, the solid lines axe for exact solutions, and the circles 
are for numerical solutions. 

Figures 1-12 are for Burgers’ equation (4.1); 

Figures 13-20 are for the linear two dimensional equation with discontinuous 
initial data (4.3), with Ax — 

Figures 21-32 are for Riemann problems for the nonconvex flux (4.4). Figures 
21-26 correspond to «i e ft = 2, u r ight = -2; Figures 27-32 correspond to «i e f t = 
—3, Vright = 3; 

Figures 33-38 are for Riemann problems for the nonconvex flux (4.5). 
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Figure 4: 4-4-LF-ENO, 
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Figure 5: 4-4-LF-ENO, Ax = , t = - i 
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Figure 8: 3-3-LF-ENO, Ax = -L- , t - 1.1 
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Figure 9: 3-3-LF-ENO, Ax = -^ 
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Figure 10: 4-4-LF-ENO, Ax 
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Figure 11: 4-4-LF-ENO, Ax = • 
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Figure 13: 3-3-LF-ENO, y = 0, t = 2 
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Figure 16: 4-4-LF-ENO, y - — 0.4, t-2 
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Figure 22: 3— 3-LF-ENO, Ax - -zrr 
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Figure 23: 


3-3-LF-ENO, Ax = -^q 
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Figure 24: 4-4-LF-ENO, Ax = 
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Figure 25: 4-4-LF-ENO, Ax = -Lr 
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Figure 26: 4-4-LF-ENO, Ax = ^ 
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Figure 28: 3-3-LF-ENO, Ax - 
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Flgure 29: 3-3-LF-ENO, Ax = 
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Figure 30: 4-4-LF-ENO, Ax = ^77 
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Figure 31: 4-4-LF-ENO, Ax = 
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Figure 35: 3-3-LF-ENO, Ax = 4 - 
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Figure 38: 3-3-EO-ENO, Ax = 
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